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Abstract 

We explain in detail how to estimate mean values and assess statis- 
tical errors for arbitrary functions of elementary observables in Monte 
Carlo simulations. The method is to estimate and sum the relevant 
autocorrelation functions, which is argued to produce more certain 
error estimates than binning techniques and hence to help toward a 
better exploitation of expensive simulations. An effective integrated 
autocorrelation time is computed which is suitable to benchmark ef- 
ficiencies of simulation algorithms with regard to specific observables 
of interest. A Matlab code is offered for download that implements 
the method. It can also combine independent runs (replica) allowing 
to judge their consistency. 
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1 Formulation of the problem 



In this article we collect some theory and formulas for the practical estima- 
tion of statistical and systematic errors in Monte Carlo simulations. Em- 
phasis is put on the estimation of in general nonlinear functions ('derived 
quantities') of the primary expectation values. The strategy focuses on the 
explicit determination of the relevant autocorrelation functions and times. 
We shall discuss why this is advantageous compared to the popular binning 
techniques, which handle autocorrelations only implicitly. A Matlab code is 
made available, which implements the method described here. 

The material given here is only partially new. It has been briefly discussed 
in the appendix of [T] and in internal notes of the ALPHA collaboration [2] 
and in lecture notes [3]. However, due to the fact that mastering this topic is 
an indispensable prerequisite for every student or researcher in the popular 
field of lattice simulation, we felt that there is a gap in the readily citeable 
literature which we want to fill here. 

We assume to have a number of primary observables enumerated by a 
Greek index with exact statistical mean values Aa. The object to be esti- 
mated is a function of these, 

F^fiA,,A2,...)^fiA^). (1) 

A typical example is given by reweighting, where a quotient F = A1/A2 has 
to be determined. A more involved case are fit-parameters determined by 
correlation functions over a range of separations. 

For the estimation we employ Monte Carlo estimates of the primary ob- 
servables given by a^. For each observable there are i = 1,2, . . . , N successive 
estimates separated by executions of some valid update procedure. We as- 
sume that the Markov chain has been equilibrated before recording data 
beginning with i = 1. This has to be secured from case to case, and we in 
particular recommend visual inspections of initial histories of observables to 
decide on safely discarding data. A posteriori one should also check that the 
equilibration time was much longer than autocorrelation times found in the 
analysis of the supposed equilibrium data. 

A key theoretical quantity for error estimation is the autocorrelation |H [5] 
function Ta/3 defined by 

<K-^«)(ar-^/3)> = ra/3H. (2) 
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that correlates the deviation of the i'th estimate for with the deviation for 
variable f3 after performing n > updates. Averages ( ■ ) here and through- 
out this paper refer to an ensemble of identical numerical experiments with 
independent random numbers and initial states — a concept useful for er- 
rors that are themselves statistical in nature. The dependence of Tap on the 
separation n alone is due to being in equilibrium. 

We take the simulation to produce configurations (f)i distributed with 
the normalized Boltzmannian P{4>)- The update algorithm is specified by 
transition probabilities W{(j) — >■ 0') and the measurements by = Oa{4>i)- 
Then the autocorrelation function in equilibrium is given as a double sum 
over configuration^, 

r«M^) = - <P'){OM) - A.){Op{<P') - Ap). (3) 

Here characterizes transitions by a succession of n update steps. If W 
obeyed detailed balance with respect to P, we would conclude symmetry 
in the indices ap. In practice we usually have only stability and no such 
symmetry seems to be implied in general. It is however natural and useful 
to define 

Tapi-n) = r^a(n) (4) 

implying a combined symmetry. Obviously rQ,^(0) is the symmetric covari- 
ance matrix with the ordinary variances on the diagonal. It is a (static) 
property of the simulated statistical system alone, while F-values for nonzero 
separation are dynamic and depend on the Monte Carlo algorithm in use. 

We end this theoretical setup by a slight generalization of our analysis. 
We add another index to the estimates of the primary observables, a^*". The 
index r = 1,2, ... ,i? counts a number of statistically independent replica. 
These are independent simulations which typically arise from repeated runs, 
or in particular, if parallel computers are used in the 'trivial parallelization' 
model. For each replicum r there are i = 1,2, Nr successive equilibrium 
estimates separated by executions of the same update procedure. Obviously 
([2]) becomes 

{{ai^' - A^){aj^' - Ap)) = 5,,,r,^(j - z). (5) 

^ We here use a discrete notation for sums over configurations. The generalization to 
continuous fields with integrations is obvious. 

^The reader is invited to simplify all formulas to i? = 1 in a first reading. 
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2 Errors and biases 

2.1 Primary observables 

We define the per replicum means 



^ 1=1 



and in terms of tliem tlie natural estimator for the primary quantities A, 

R 



hT^^rK (7) 

with the total number of estimates 

N = J2Nr. (8) 



r=l 



R 



r=l 



Introducing the deviations 

K = ^a- ^a, 6a = aa - Aa- (9) 

we state that these estimators are unbiased, 

(0 = = (10) 

We assume a normal distributioijfl which is then completely defined by the 
covariance matrix given by 

1 1 

( ^« ^1 ) = ]^ E ^"/^(J - ^)Srs = J^Cap 5rs X (1 + 0(r/iV,)) (11) 
i,j=l 

where 

oo 

Cap = J2 ra/3(t) (12) 
t=—oo 



^ Note that the a^f themselves need not be Gaussian distributed and that large Nr 
makes the Gaussian by the central limit theorem. 
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does not depend on the run-length. Here and in the following it is assumed 
that a finite scale r characterizes the asymptotic exponential decay of F^^ 
with |t|. The claimed A^^- dependence holds for runs with lengths Nr ^ t 
and if this condition is violated, an error estimation is hardly possible. For 
the covariance matrix of 

{Lh) = ^Cap X (1 + 0{Rr/N)) (13) 

follows. Roughly speaking, as expected, aa differs from by an error of or- 
der 1 / ^/N and the task of an error analysis amounts to a reasonably accurate 
estimation of Cq,^ which includes all autocorrelation effects. 



2.2 Derived quantities 

To determine F we consider estimators 

F = /(a.) (14) 

and 

1 ^ 

^^N^Z^rfiK). (15) 
r=l 

We assume the estimates of the primary obscrvables to be accurate enough 
to justify Taylor expansions of / in the fluctuations, for example 



F^F + J2faSa + lj2j'o.pSjf3 + ... (16) 
a 

with derivatives 



2 

a/3 



df d'^f 
" OA.' -^"^ " dA^dAfs ^^^^ 
taken at the exact values Ai,A2,.... It follows that our estimators for F are 
biased unless / is linear, 

{F-F)^^Y.f'^^^'^P^^i^-^)^ (18) 

Usually, this bias is negligible compared to statistical errors due to large 
enough N. Replica (i? > 2) can however be used to control and, if deemed 
to be appropriate, cancel the leading bias by using 

{P-F)c^^{F-F) (19) 
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to replace 

The error ap is to leading order given by 

a].= {{F-Fr)c^^Cp (21) 

with 

Cp = ^ faff3 Ca/3 ■ (22) 
a/3 

We rewrite this expression as 

(Tp = vp (23) 

with the effective 'naive' - i.e. disregarding autocorrelations - variance rele- 
vant for F 

fF = $^/a//3r„^(0) (24) 

a/3 

and the integrated autocorrelation time for F 

^ oo 

W = ^ E J^f'^fMt)- (25) 

-'^ t=-oo a/3 

This one number encodes the efficiency of the algorithm in use for a deter- 
mination of F, if its execution time per update is known. To appreciate 
this notation for the ratio of the true to the naive error two extremal cases 
are instructive. In the absence of autocorrelations, T a/sit) oc St ft, we have 
2Tint,F = 1 and fl23l) is the standard error. For a purely exponential be- 
haviour, Taisit) OC exp(— |t|/r), this scale is recovered, rint,F = t + 0(r~^). 
In general there are however many contributions and rint,F is more like a 
'dynamical susceptibility'. Another interpretation of fl25]) is, that only the 
reduced number of N/ (2rint,F) estimates are effectively independent to bring 
down statistical errors. 

If we have R>2 replica another consistency check is useful. The replica 
estimates /(a^) are assumed to be normally distributed. Their variances are 

{ifiK)-F)') = ^. (26) 
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We may imagine to do a fit of the replica estimates to a constant K by 
minimizing 

r ' 

The minimum occurs precisely aX K = F given in (fT5|) . The probability to 
find in such a 'fit' a minimal > 2; is given by 

g = l-P(x/2,(i?-l)/2), (28) 

where P is the incomplete Gamma function, 



1 



u 



pmn) = -^l duu""-^ exp(-u). (29) 

r(n) Jo 

Taking the actually observed x = x^i^) we obtain the so-called goodness of 
fit. It is very unlikely to see Q- values much less than 0.1, if the underlying 
distribution has the required properties, and this is hence a good consistency 
check In a histogram plot we may in addition inspect 

P = ^^'"^ ~ ^ (30) 
" ap^N/Nr - 1 

which should have zero mean and unit variance. 

In the next section we develop numerical estimators for vp and Tjnt.F, 
which at present are not yet of practical use, since they so far involve the 
numerically unknown exact Aa and Tap- 



3 Extraction of errors from measured data 
3.1 Error estimators 

We now describe what we would like to call the F-method of error estima- 
tion (as opposed to binning, see App|B|, where we explicitly estimate the 
autocorrelation function. 

In terms of the given estimates a^^' we define an estimator Tap of the 
autocorrelation function 

R Nr~t 

= -J^TTWi E E « - - (31) 
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where < t <^ Nr for all r is understood. An important principle underlying 
this choice is to avoid unnecessary noise from terms which vanish (only) 
on average. Therefore the correlation is first taken within each replicum and 
only then averaged. Using arguments of the same type as for ffTTl) the leading 
bias of this estimator is evaluated, 

(f.,(t))-r.,(t)^-^. (32) 

It comes from subtracting in (!3T!) the ensemble mean instead of the exact 
average. For uncorrelated data this is the well known bias in the estimator 
of the variance cancelled by the — 1 instead of normalization. For the 
moment we neglect this bias of order 1/A^ relative to the l/\fN statistical 
error in Tap to be worked out in the next subsection, and come back to it 
later. 

For the derived quantity F only the projected autocorrelation is required, 

fF(t) = 5^/JA/3W (33) 



where /„ is the gradient as in flTTl) but now evaluated at arguments ai, 02, . . ., 
which inflicts a small relative error of order 1/A^ that we do not attempt to 
cancel. 

As practical estimators for vp and Cp we take 

= fi.(0) (34) 

and 



Cf{W) 

t=i 

A bias results from the truncation of the autocorrelation sum, 



w 



fp(0) + 2^fi.(t) 



(35) 



i^MEl)^._e.p(-^/.). (36) 

For a purely exponential Tp this is an approximate equality, otherwise an 
order of magnitude statement with the scale r introduced at the end of 
section 12. 1[ The summation window W has to be chosen with care. On 
the one hand it has to be large compared to the decay time r for a small 
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systematic error, on the other hand with too large W one includes terms with 
negligible signal but excessive noise. In the next subsection we shall define 
an automatic procedure to choose W. 

In a frequently met case we have many primary observables but only 
want to estimate a few functions of them. Then it is advantageous to actually 
not compute the full autocorrelation matrices Ta/sit) but to just form and 
fa and immediately project the measured data 

< - «r = E ^< (37) 

a 

to base the remaining autocorrelation analysis directly on them. 

Another practical point concerns the computation of fa- It may well be 
inconvenient or even impossible to analytically specify the gradient of / and 
we may prefer to evaluate it numerically based on the supplied routine for / 
itself only. A natural scale ha is extracted from the data, 

ha = ^J^^ , (38) 

and we take 

la ~ ^[/(^i' a2,...,aa + ha,...)- /(oi, Os, . . . , a„ - . . . )] (39) 

as a numerical estimate for the gradient of / with negligible errors of 0{h'^) oc 
1/N. 

3.2 Error of the error 

When an error is assessed numerically in a Monte Carlo it is itself affected 
by statistical error. This is usually not quantified. But to rate algorithms on 
the basis of the Ti^t,F they produce for physically interesting F it is better to 
at least have a rough idea. This is even more true, if one wants to falsify a 
theoretical prediction based on one of these '2.5 cr discrepancies', while maybe 
the estimate of ap itself is uncertain by a factor 2. 

In appendix [X] the very convenient approximate formula 

{{CAW)-C,r)^^-^^^^Cy (40) 
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is derived that was to our knowledge first given in For the quotient 
we similarly fincfl with the help of §3) and ([58]) 



{{fintAW) - Tint.F)') + nnt,F)rL,F- (42) 

Beside the statistical error of Cf{W) there is the systematic error from 
the ly-truncation in fl36l) . We define as optimal the W value that minimizes 
the sum of the absolute values of these errors and yields a total relative error 





for the final error estimate 



- min I exp {-W/t) + 2^W/n] (43) 
2 w \ / 



= (44) 

evaluated at the optimal W. It is a function of N/ r and plotted in Fig{T] for 
relevant arguments. The minimum leads to a simple transcendental equation 
that we solved numerically, but it can also be well approximated hj W = 
rln(iV/r)/2. The asymptotic behaviour is 6tot oc ^ln{N)/N. Due to the 
exponential behaviour fl36|) the systematic component of the error becomes 
negligible at large N, albeit only slowly, 



4ys(0"F) T 



4tat(^F) 2W \n{N/T) 



(45) 



A simple and popular alternative to the analysis method just described 
are binning methods, where data are pre-averaged over sections of the history 
of length B. These bin averages are taken as independent estimates. For both 
ordinary and jackknife binning (see appendix [B] for more details) a relative 
systematic error of the error estimate occurs, which decays only proportional 
to t/B. The relative statistical error of the error due to the finite number 



^ We here use error propagation and the upper bound in ([55]) . This leads to an upper 
bound on the error of fint.F regarding the terms beyond the one oc W which is both exact 
and dominant. 
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Figure 1: Relative total error of the error estimate for the F-method and for 
binning. 

of bins is yJlB jN . By balancing these errors similarly as before we find the 
total error at optimal bin size 



O'Fbin 2 B 

which is reached for 

B = r(2iV/r)i/3 (^47) 

and is also shown in FigUl Here the ratio of systematic to statistical error is 
constant, 

(48) 



^sys(o'Fbm) _ 1 



4tat(<7Fbm) 2 
11 



In comparing the r-method with binning it is clear now that the bin size 
B and the window size 2W play a very similar role. All advantages of the 
F-method derive from the exponential in W rather than inverse power in B 
behaviour of the systematic error of the error. One could say, that in this 
way we have an improved estimator of the error whose own error decays with 
a (slightly) larger power of A^. 

We briefly come back to the bias in our error estimate caused by the right 
hand side of (l32l) . Cancelling this term when estimating Cp by (135|) amounts 
to enlarging 

/ 2W + 1\ 
CAW) - C,{W) i^l + , (49) 

which at the optimal W value represents a correction down by a factor 1 / -\/iV 
(up to logs) relative to both contributions in fHSj) . Since the correction is 
simple and known in closed form we include it in applications. 



3.3 Automatic windowing procedure 

We present now an algorithm to automatically choose W that is similar but 
not identical to the one proposed in [1]. To determine W in an actual data 
analysis, we start from a hypothesis that r ~ STi^t,F with some factor S. 
More precisely, we solve 

oo 

^fintAw) = (^M-s\t\/f{w)) (50) 

t= — OD 

for r, to get a formula that is more precise for small f , 

s _ / 2fi„t,F(w^) + n _ 1 1 , . 

If T^int F < 1/2, we set f{W) to a tiny positive value. Then we compute for 

= 1,2,... 

g{W) = exp[-W/f{W)] - f{W)/VWN . (52) 

Up to a factor this is the VT-derivative defining the minimum in fH3l) . The 
first value W where g{W) is negative (sign change) is taken as our summation 
window which under our hypothesis is self-consistently close to optimal. For 
most practical systems, algorithms and observables r and rint.F are of the 
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same order, and S" = 1 ... 2 is a reasonable choice. A too large S leads to (in 
principle) unnecessarily large statistical errors of the error. 

It is important to verify by eye, at least for a representative set of ob- 
servables, that fint,F(W) with its statistical errors fH2|) as a function of W 
exhibits a plateau around the automatically chosen value. If this is not sat- 
isfactory, one has to modify S to achieve this. A similar procedure could be 
set up to semi-objectively choose the bin size B although this will be less 
stable, since the plateau behaviour will be less pronounced. 

4 Conclusions 

We have given a detailed description of the error analysis of Monte Carlo 
data. We found advantages in explicitly analyzing autocorrelation functions 
rather than employing binning strategies. They are due to a component in 
the systematic error of the error estimate from binning which is proportional 
to t/B. It is a finite bin-size effect in the integrated autocorrelation time 
from these procedures, which is avoided in the F method. The discussion 
includes a careful analysis of the errors of the error estimate, both statistical 
and systematic, and a minimization of their sum. Their approximate size 
can be read off from Fig{T] once one has an at least rough estimate of the 
autocorrelation decay scale r. 

A Matlab routine that implements the whole method, including the nec- 
essary plots, is described and offered for download in the internet. It is also 
capable of handling multiple runs (or subdividing one run into several parts) 
and to judge the compatibility of the multiple estimates relative to the over- 
all error estimate. The routine adapts itself mostly automatically to the 
autocorrelation times involved. This works under reasonable assumptions 
mentioned in the text, but an additional visual inspection of the autocorre- 
lations for at least part of the observables is always recommended. 

Acknowledgements. I would like to thank the members of ALPHA 
for discussions and numerous tests of the software described here and its 
precursors. Philippe de Forcrand, Andreas Jiittner, Bjorn Leder, Francesco 
Knechtli, Thomasz Korzec, Martin Liischer, Juri Rolf and Rainer Sommer 
helped with feedback on the manuscript. 
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A The Madras- Sokal formula for the statisti- 
cal error of the error 

To prove fj40l) we need to investigate 

(f„^(s)f^5(t)) = (53) 

^ N-s ^ N-t 

^ i=i j=i 

This formula is for only one replicum at the moment to avoid even more 
indices. The essential approximation consists now of approximating this 
four-point function by its disconnected part, 

( f^^{s)f,s{t) ) ^ {T^p{s) - Co,p/N){T,s{t) - C,s/N) + (54) 
1 ^ 

+]V2 5Z " ^)r/35(j -i + t-s) + TasU - i + t)Tp^{j - i-s) 

This is plausible for large A^. At least gathering from experience with or- 
dinary statistical systems the connected part makes a smaller contribution 
since all arguments have to be close to each other. We have also assumed 
\s\,\t\ <C to extend the sums with relative errors of order |s|/A^, \t\/N. 
Furthermore we use for a rapidly decaying function g the approximation 

i,j=l k=—oo 

similarly as in (fT2!) to obtain 

^ oo 

((f^(s)-r^(s))(f^(t)-r,.(t))) = - J2 ^Fm^F{k+t-s)+TF{k-t~s)]. 

k=—oo 

(55) 

Summing now t,s = —W, . . . ,W and assuming that each factor Tp cuts off 
the sum unless its argument is much smaller in absolute value than W we 
derive 

2 °° 2 

{(CpiW)-Cpf)^-i2W + l) J2 rpik)Tpil) = -i2W+l)Cl (56) 

fe,/=— oo 
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with corrections of 0{t/W). 

By very similar steps one may also derive 

( - v^){dFiW) -C^))^ ^Cl. (57) 
A more complicated case is 

9 °° 2 
( {v, - v,f )^j^Y. (r^fc))' < ^VfC,, (58) 

fc=— oo 

where the inequality certainly holds if < r(A;) < r(0) is true. Even if this 
happens to be violated out in the tail of F, we still expect the inequality to 
hold for the sum (dominated by small \k\) in practical applications. 

Repeating the above analysis with several replica leads to the same for- 
mulas and thus yields (HO]) . In view of the approximations we note that the 
practical experience is, that the resulting errors of Tint,F are consistent with 
their scatter when several runs are made. 



B Binning methods 

Here we give the basic formulas for error estimation by binning. We assume 
data a\ where possible replica are sewed together to one history of length 
= BNb which we divide into Nb sections of B consecutive measurements 
each0. We form bins or blocked measurements 

hi = ^Y^a^t'^^^\k = l,...,Ns (59) 
j=i 

The mean values are all identical h = a = a. The idea is now to choose 
B large enough such that the bins 6^ can be regarded as approximately 
uncorrelated, but small enough such that there remains a large number Nb 
of such 'events'. As an error estimator one then takes 

Nb 
k=l 

^ The breaks of autocorrelation at replica boundaries are neglected and a possible 
remainder is discarded. 
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If we Taylor-expand / again, the relevant correlation matrix is 



(61) 



u,v 



^(r«^(M -i) + r„^(i -u)) + ^Y^ r„^(i - j) 



u,i i,j 



where indices u,v run through bin k only while i,j range from 1 to as 
before. The result is dominated by the first term 



while the other terms approximately sum to —Cap/N. The coefficient oir/ B 
is exact if Tap is simply exponential, but in practice this should be read as 
0{t/B). This term is the analog of (136|1 for the F- method. We see that 
the integrated autocorrelation time is effectively constructed by a double 
sum within one bin of size B and the t/B error is a finite size effect which 
masks the actual exp(— 5/r) expected from truncation. The many bins are 
only used to tame statistical fluctuations of the estimate but do not improve 
the systematic error. The analogous flnite size effect for the F-method occurs 
only at order t/N. A cancellation of the t/B or exp(— H^/r) effects dominant 
for the two methods was not found to be practical, as one would have to use 
a noisy estimate for r obtained from the data and would not know the precise 
pre-factor. 

In ( 160|) the function / has to be evaluated on averages over only B events, 
that is only the Ns^th fraction of all measurements. This may be affected 
with stability problems in the case of flts for example. Compared to our 
earlier discussion, the gradient of / is here implicitly formed stochastically 
with larger average spreads ~ ^/Nb ha- Both problems are cleverly overcome 
by going to complementary or jackknife bins. 



Lj2Tap{u~v)c^^Cap{l-T/B) 



(62) 



u,v 




(63) 



An elementary calculation yields the relation 



( 




)(cJ-«/3)> 



1 



2 ((&^-a,)(6j-a^)>. 



(64) 



'a 



{Nb - 1) 
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Hence the resulting jackknife error estimator is 



AT — 1 

^ k=l 

Here all evaluations of / practically work on the full statistics, and the spread 
of its arguments is of the same order ha as taken for the numerical derivative 
in the F-method. The result of the above analysis of systematic errors is 
however unchanged for the jackknife estimate. 

For completeness we finally mention that also under binning the scatter 
between F, the function / applied to the mean over all data, and the average 
of /(c^) over all bins can be used to cancel the leading 0(1/A^) bias in the 
estimate of the mean of derived observables. In analogy to fl20l) we have to 
substitute 

. Nb 

F-.NsF+{l-Ns) — J2fi<^^)- (66) 

^ fe=i 

C An implementation of the F-method 

In the ALPHA collaboration we have recently performed most data analyses 
in Matlab, since we found that it combines comfortable programming, a 
robust library, interactive graphics and acceptable speed suitable for this 
purpose. An improved version of the well tested core-routine is described 
in the first subsection. While it has been applied to many realistic data 
sets in the ALPHA projects, we here develop a simulator, that allows for 
clinical testing with exactly known errors in a situation that we consider 
representative for simulations in lattice QCD. 

C.l Description of the Matlab routine UWerr 

The calling sequence of the function is 

[value, dvalue, ddvalue, t auint jdtauint ,Q] = ... 
UWerr(Data,Stau,Nrep, Name, Quantity, PI, P2, . . .) 

with the meaning of the input and output arguments described below. Sim- 
ilar information is obtained from help UWerr under Matlab. 

Input: 
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Data is a matrix filled with the primary estimates a^^ from R replica with 
Ni, N2, . . . , Nfi measurements each, 







1,1 

^2 


1,1 

ag 




2,1 


2,1 


2,1 
«3 


Data = 


Ni,l 


Ni,l 
«2 


Afi,l 
«3 




1,2 
a{ 


1,2 
«2 


1,2 
«3 



... \ 



Affl.fl Affl.i? Nii,R 
\ "1 "2 "3 ■ ■ ■ / 

All input arguments that follow can be omitted or set to empty (=[]) and 
are then assigned default values. 

If data after reading them into Matlab occur in a different matrix format, 
e.g. one long vector, the Matlab commands reshape may be employed to 
transform them very efficiently and safely. 

Stau is the estimate S = r/rjnt explained in section (13. 3p . The default is 
S = 1.5. It is tried to determine an optimal window W as described, but it 
is never taken larger than min{A^r}/2. If the windowing condition cannot be 
met, a warning is printed. 

Nrep = [A^^i, N2, . . . , N^i] is a vector whose elements specify the lengths of the 
individual replica. As default all data are interpreted as one replicum, Nrep 
= [N] . Even if only one run has been made, one may specify Nrep such as 
to analyze sections of the history independently and investigate the scatter 
of sub-averages. 

Name can be a string (default: 'NoName') which is used as a title in the 
generated plots (see below). If the argument is given but is not a string (an 
integer for example) no plots are generated. 

Quantity can either be an integer from the range of the index a, in which 
case the corresponding primary quantity is analyzed. Alternatively it may 
be a function handles @fun if fun is the name of the function /(Ai, yl2, .. .) 
which has to be user-defined and has as first argument a row-vector which 
has as many entries as Data has columns. In this case the corresponding 
derived quantity is analyzed. Optional further arguments PI, P2, . . . are 

^ At present also a string with the name can be entered, which however is obsolete and 
announced not to be supported in future versions of Matlab. 
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passed on to the function as second, third, . . . argument. The default value 
is Quant ity=l. 



Output: 



value, dvalue,ddvalue,tauint,dtauint,Q have their obvious meaning as 
best estimates of F,aF as defined in (|Ti).( l2Ti) . statistical error of the error 
from (l40l) . rint,F and (5(rint,F) given in (125|) and extracted from fH2l) . The 
Q-value (l28l) follows if i? > 2. Note that Cf,vf and the chosen W can be 
reconstructed from the output if necessary. 

The leading bias is corrected according to (!20|) for R > 2. If this correc- 
tion amounts to more than crir/4 the user is warned. The routine generates 
— unless suppressed — plots of our estimate for p{t) = Tp(t)/Tp{0) and 
^int,F(W^) with errorbars in the relevant range to inspect the required plateau 
behaviour. Errors of p{t) have been added in version 6 of the UWerr. They 
are computed with formula (E.IO) in [7] taking A = t + W. If there are two or 
more replica their distribution fISPl) is shown as a histogram with its Q-value 
f l25]) in the title. In addition a histogram of all a^'" is displayed in case that 
a primary observable is analyzed. 

The code (m-file) of the Matlab routine UWerr . m described here together 
with a sample function effmass.m used in the following test can be down- 
loaded from www-com.physik.hu-berlin.de/ALPHAsoft 

C.2 Test in a simulator 

As a test case we think of determining the 'effective' mass 



from a correlation G{z) = exp{—mz) in euclidean time z, which we take 
simply exponential since we are not interested here in systematic errors of 
m. To fake noisy (auto) correlated data for G we employ normally distributed 
random numbers ?7i, i = 1,2,... with zero mean and unit variance to produce 
a sequence {z/j} by recursively setting 



m = ln{G{z)/G{z + 1)) = f{G{z), G{z + I)) = F 



(67) 



^1 = Vi, 



a2 r]i+i + avi 




with 



2r- 1 



(69) 



a = 



2r + 1 



> -. 

- 2 



^ I thank Hubert Simma for pointing out an improvement here. 
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normalized autocorrelation of mass 




-0.5 



5 



X. , with statistical errors of mass 




Figure 2: Generated plots of TpiJ:) and Tiot,F{W) with a vertical line at the 
ly- value picked automatically and a horizontal one at the resulting estimate 

of Tint.F- 



A short calculation yields 

oo 

( i^i ) = 0, ( i^i+ti^i ) = a\ al*l = 2r. 



(70) 



t=—oo 



We use three independent such sequences , vf'\ vf' and choose for them 
integrated autocorrelation times r^^^ < r^"^^ = t^^\ 

Synthetic data to measure m = f{Ai, A2) are formed, 



al = G(0)+g(z.f) + z.f), a^, = G{1) + qiu^ + u^') ■ (71) 



(1) , ,,(3)x 
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replica distribution (mean=0,var=1) for mass » Q = 0.28 « 



1.5 



0.5 






_2 -1.5 -1 -0.5 0.5 1 1.5 2 



Figure 3: Replica distribution of {pr} with the Q- value in the title. 



The real parameter q sets the strength of the noise, and the presence of u^^^ 
in both observables creates a correlation between them, as in real life. 

All relevant quantities of the error analysis can be worked out exactly for 
this case. 



VAi = VA2 = 2g , Tint.Ai = 1"int,A2 = (t^ ' + t'^ ') /2 

vp = 2g^(l + exp(2m) — exp(m)). 



Tint,F 



^(1) + ^V2 + l ^(2) 



m? + 1 



m? + 1 



(72) 
(73) 
(74) 



with 



m = 2 sinh(m/2). (75) 

Note that in the weighed sum giving Tint,F the time scale r*-^-* contributes 
with a small weight, due to correlations between a\,a^ and the cancellation 
of the corresponding fluctuations in the effective mass. 

For our test we now set m = 0.2, r^^^ = 4, r^^^ = 8,q = 0.2 which leads to 
vp = 0.1016, Tint.F = 7.92. We generate 8 x 1000 estimates corresponding to 
an error in the mass of 0.0142 and get 
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Figure 4: Estimation of Ti^t,F by jackknife binning. A vertical line shows the 
optimal bin size (H71) . 

» Nr=ones(8,l)*(N/8) ; 

» [value, dvalue,ddvalue,tauint,dtauint,Q]= ... 

UWerr(Data, 1 ,Nr , 'mass' ,@ef f mass , 1 , 2) ; 
» [value , dvalue , ddvalue , tauint , dtauint , Q] 
ans = 

0.2128 0.0134 0.0008 7.2909 0.7396 0.2827 

The generated plots are shown in Figl2] and in Figj3l By repeating such 
an analysis 20000 times we knock down the statistical error of the error far 
enough to exhibit a systematic error of ~ —0.5% which is of the expected 
order |exp(— l^/r) with W = 37 and r = 8, while, of course, the actual 
automatic W and estimated r-values fluctuate. 

The jackknife analog of the lower part of Figj2]can be seen in Fig JH for the 
same data. It looks consistent with a less pronounced plateau, as expected. 

More relevant are tests on real data. The routine passed a number of 
such tests, but, as usual with software, with probability one some future 
application will force further evolution. 
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